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Abstract 

In the mathematical modelling of compactional flow in porous media, the constitutive relation 
is typically modelled in terms of a nonlinear relationship between effective pressure and porosity, 
and compaction is essentially poroelastic. However, at depths deeper than 1 km where pressure is 
high, compaction becomes more akin to a viscous one. Two mathematical models of compaction in 
porous media are formulated and the noninear equations are then solved numerically. The essential 
features of numerical profiles of poroelastic and viscous compaction are thus compared with asymp- 
totic solutions. Two distinguished styles of density-driven compaction in fast and slow compacting 
sediments are analysed and shown in this paper. 
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1 INTRODUCTION 

Density-driven compaction in porous media such as sediments is an important process, 
which may occur in sedimentary basins where hydrocarbons and oil are primarily formed. 
The modelling of such density-driven flow is thus important in the oil industry as well as 
in civil engineering. One particular problem which affects drilling process is the occasional 
occurrence of abnormally high pore fluid pressures, which, if encountered suddenly, can 
cause drill hole collapse and consequent failure of the drilling operation. Therefore, an 
industrially important objective is to predict overpressuring before drilling and to identify 
its precursors during drilling. An essential step to achieve such objectives is the scientific 
understanding of their mechanisms and the evolutionary history of post-depositional 
sediments such as shales. 

Fine-grained sediments such as shales and sandstones are considered to be the source 
rocks for much petroleum found in sandstones and carbonates. At deposition, sediments 
such as shales and sands typically have porosities of order 0.5 or 50%. When sediments 
are drilled at a depth, say 5000 m, porosities are typically 0.05 ~ 0.2 (5% ~ 20%) [1]. 
Thus an enormous amount of water has escaped from the sediments during their deposi- 
tion and later evolution. Because of the fluid escape, the grain-to-grain contact pressure 
must increase to support the overlying sediment weight. Dynamical fluid escape depends 
lithologically on the permeability behavior of the evolving sediments. As fluid escape pro- 
ceeds, porosity decreases, so permeability becomes smaller, leading to an ever-increasing 
delay in extracting the residual fluids. The addition of more overburden sediments is 
then compensated for by an increase of excess pressure in the retained fluids. Thus over- 
pressure develops from such a non-equilibrium compaction environment [2]. A rapidly 
accumulating basin is unable to expel pore fluids sufficiently rapidly due to the weight of 
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overburden rock. The development of overpressuring retards compaction, resulting in a 
higher porosity, a higher permeability and a higher thermal conductivity than are normal 
for a given depth, which changes the structural and stratigraphic shaping of sedimentary 
units and provides a potential for hydrocarbon migration. 

Compaction is the process of volume reduction via pore-water expulsion within sedi- 
ments due to the increasing weight of overburden load. The requirement of its occurrence 
is not only the application of an overburden load but also the expulsion of pore water. 
The extent of compaction is strongly influenced by burial history and the lithology of 
sediments. The freshly deposited loosely packed sediments tend to evolve, like an open 
system, towards a closely packed grain framework during the initial stages of burial com- 
paction and this is accomplished by the processes of grain slippage, rotation, bending and 
brittle fracturing. Such reorientation processes are collectively referred to as mechanical 
compaction, which generally takes place in the first 1-2 km of burial. After this ini- 
tial porosity loss, further porosity reduction is accomplished by the process of chemical 
compaction such as pressure solution at grain contacts. It is worth pointing out that 
consolidation is a term often used in geotechnical engineering and implies the reduction 
of pore space by mechanical loading. The fundamental understanding of mechanical and 
physico-chemical properties of these rocks in the earth's crust has important applications 
in petrology, sedimentology, soil mechanics, oil and gas engineering and other geophysical 
research areas. In spite of its geological importance, the mechanism leading to pressure 
solution is still poorly understood [3]. 

The main aims in this paper are to determine and compare the essential features of 
the poroelastic and viscous compaction in a comprehensive way and to understand these 
mechanisms by using new asymptotic solutions and the comparison with full numerical 
simulations as well, which will greatly extend the earlier work [2-4]. Another primary 
concern of this paper is to try to formulate a new and more realistic visco-poroelastic 
compaction relation. 



2 MATHEMATICAL MODEL 



For the convenience of investigating the effect of compaction in porous media due to pure 
density differences, we will assume the basic model of compaction is rather analogous to 
the process of soil consolidation. The porous media act as a compressible porous matrix, 
so that mass conservation of pore fluid together with Darcy's law leads to the 1-D model 
equations of the general type [3,4]. Let t be time and z be the space co-ordinate directing 
upwards, the governing equations can be written as 

d[ps(1 m (f))] +^- z [p s (l-<P)u s ] = 0, (solid phase) (1) 

^ + ^ = 0, (liquid phase) (2) 

4>{u l - u s ) = ^y-[G^ - ( Ps - Pl )(l - fig], (Darcy'slaw) (3) 

where is the porosity of the pores saturated with water. u l and u s are the velocities 
of fluid and solid matrix, k and p are the matrix permeability and the liquid viscosity, 
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pi and p s are the densities of fluid and solid matrix, p e is the effective pressure, G is a 
constant of the properties in porous media, and g is the gravitational acceleration. In 
addition, a compaction relation is needed to complete this model [4,5]. By assuming the 
densities p s and pi are constants, we can see that only the density difference p s — p\ is 
important to the flow evolution. Thus, the compactional flow is essentially density-driven 
flow in porous media. 



2.1 Poroelasticity and Viscous Compaction 

Compaction relation is a relationship between effective pressure p e and strain rate e = ^ 
or porosity <p [6]. The common approach in soil mechanics and sediment compaction is 
to model this generally nonlinear behaviour as poroelastic, that is to say, a relationship 
of Athy's law type p e = p e (0), which is derived from fitting the real data of sediments. 
Athy's poroelasticity law is also a simplified form of Critical State Theory. A common 
relation representing the poroelasticity is 

Dt s dz' Dt dt + dz 1 {> 

where K s is a modulus of sediment compression. As p s is a constant and can thus be 
eliminated by multiplying equation (1) by l/p s , and we get 



d{i-<t>) uS d{i-<t>) or i sq-fl du* 

dt +U dz ° r 1-0 Dt ~ dz> (5j 

combining with the previous equation (4), we have 

Pe=Pe(<j)), (6) 

which is the Athy's law for poroelasticity. However, this poroelastic compaction law is 
only valid for the compaction in porous media in the upper and shallow region, where 
compaction occurs due to the pure mechanical movements such as grain sliding and 
packing rearrangement. In the more deeper region, mechanical compaction is gradu- 
ally replaced by the chemical compaction due to stress-enhanced flow along the grain 
boundary from the grain contact areas to the free pore, where pressure is essentially pore 
pressure. A typical process of such chemical compaction in sediment is pressure solution 
whose rheological behavior is usually viscous, so that it sometimes called viscous pressure 
solution or viscous creep. 

The mathematical formulation for viscous compaction is to derive a relation between 
creep rate e and effective stress a e . Rutter's creep relation is widely used [7,8] 

A k c wD gb 

6 = -p-d^^ (7) 

where a e is the effective normal stress across the grain contacts, A k is a constant, Co is the 
equilibrium concentration (of quartz) in pore fluid, p, d are the density and (averaged) 
grain diameter (of quartz). D g b is the diffusivity of the solute in water along grain 
boundaries with a thickness w. Note that a e = —Gp e and e = With this, (7) 

becomes the following compaction law 

p e = -£V.u°, g= Gpsd3 . (8) 
A k c wD gh 
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More generally speaking, £ is also a function of porosity 0. The compaction law is 
analogous to the viscous compaction laws used in studies of magma transport in the 
Earth's mantle [9,10]. 

2.2 Boundary conditions 

The boundary conditions for the governing equations are as follows. The bottom bound- 
ary at z = is assumed to be impermeable 

u s = u l = 0, (9) 
and a top condition at z = h is kinetic 

h = rh s + u s , (10) 
where m s is the sedimentation rate at z — h. Also at z = h, 

= 00, Pe=P0, (11) 

where po is the applied effective pressure at the top of the porous media, and 0o is the 
initial porosity. 



3 Non-dimensionalization 

If a length-scale d is a typical length [8] defined by 



and the effective pressure is scaled in the following way 

= G{p e - po) 
(p s - Pi)gd' 

so that p = 0(1). Meanwhile, we scale z with d, u s with rh s , time t with d/m s , perme- 
ability k with k . By writing k(<f>) = k k*, z = dz*, and dropping the asterisks, we 
thus have 

-^ + l:[(WK] = o, (14) 



dt dz 
dt dz 



d<P + ^ = 0, (15) 



A M0)h£-(l-0)]. (16) 



[ dz 

The poroelastic relation becomes 

p = p((f>) (17) 

and the viscous relation is 

du s , , 

P=-^. (18) 
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where 

A = k ^~P^ . (19) 
lim s 

Adding (14) and (15) together and integrating from the bottom, we have 

u s = -<j)(u l -u s ) = -u, (20) 
where u = 4>{u l — u s ) is the Darcy flow velocity. Now we have 

£ + |[<i -*)«] = ». <«) 

« = -«W[^-(l-*)l- (22) 
The constitutive relation for permeability k((j>) is nonlinear [11], and its typical form is 

W) = (-jr)™ rn = 8. (23) 

Different formulations of compaction relation may lead to different compaction models. 
One way is to use a relationship between effective pressure p and matrix velocity u s (or 
u s in 3-D form) as given in (8). However, a more common way is to write a relation 
between p and porosity 0. Formulating the compaction relation in this way, we have 
Poroelastic Model: 

-<ng -a-«]>. (-) 

p= -[In ^-(0 o -0)], (25) 
a 

which is a relation of Athy-type. a = 0(1) is usually called the compaction or consoli- 
dation coefficient. The boundary conditions are 

j|-(l-0) = O, at z = 0, (26) 
= 0o, h = m(t) + \A m &-(l-<j>)] at z = h(t). (27) 

00 02 



Viscous Model: 



"=4«s )m[ i- (1 - m (29) 



The boundary conditions are 



^ - (1 - 0) = 0, at * = 0, (30) 
>o, h = rh(t) + A A m [^ - (1 - 0)] at z = Zi(t). (31) 

00 
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where m(t) = 0(1) is a prescribed function of time, which can be taken to be one for 
constant sedimentation on top of the porous media. Obviously, m — if there is no 
further sedimentation and no increasing loading on top of the porous media. 

It is useful for the understanding of the solutions to get an estimate for A by using 
values taken from observations and earlier work [1, 4, 11]. By using the typical values 
of pi ~ 10 3 kg nT 3 , p s ~ 2.5 x 10 3 kg m" 3 , k ~ 1(T 15 — 1(T 20 m 2 , p ~ 1(T 3 Nsm 2 ,(~ 
1 x 10 21 N s m" 2 , rh s ~ 300 m Ma -1 = 1 x 10 _11 m s -1 , g « 10m s~ 2 , G « 1; then 

A 0.01 1000 and d ~ 1000 m. Therefore, A = 1 defines a transition between the 

slow compaction (A << 1) and fast compaction (A >> 1). The parameter A , which is 
the ratio between the permeability and the sedimentation rate, governs the evolution of 
the pore pressure and porosity in sedimentary basins. High sedimentation rate may gives 
rise to excess pressures even in the basins with moderate permeability. 



4 Numerical Simulations and Asymptotic Analysis 
4.1 Numerical Method 

In order to solve the highly coupled non-linear equations, an implicit numerical difference 
method is used [12]. Substituting the expression for effective pressure p into the 
equation, the essential equation for porosity becomes the standard non-linear parabolic 
form 

<lH = F{z,t,<j>)<j> xz + g{z,t,<l>,<j> z ). (32) 
The first stage gives n+1 / 2 as a solution of the following equation 

+g ( Zt X +1/2 ,<t>?,-^5 z tf), (33) 

where <5 2 0j = (<f>i+i — 20j + 4>i-i) an d <5 2 0i = (l/2)(0 i+ i — At and Az are the time 

and space increments after discretisation, respectively. The second stage gives as a 
solution of the following equation 

= (^^)^(-^ n+1/2 ) 0r 1/2 )^ 2 (0r +1 +c) 

+g(z l X +1/2 :] rz S z <pr 1/2 ). (34) 

The convergence is second-order in space for this method, and 0(At) 2 ~ e in time, where 
e is a small number less than 1/2. 

The computational convergence of the calculation of this method has been tested by 
1) changing the number of grid per unit (1/Az) from 5 to 1000 in space and I /At from 
10 to 5000 in time, and by 2) comparing with the results of asymptotic results. The 
changes of grid intervals all result in the same converged results which conform well to 
the asymptotic solutions. This shows that this method is robust for the solution of the 
equations encountered in our problems. 
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4.2 Numerical Results 

We used a normalized grid by employing the rescaled height variable Z = z/h(t) in a 
fixed domain, which will make it easy to compare the results of different times with 
different values of dimensionless parameters in a fixed frame. This transformation maps 
the basement of the basin to Z = and the basin top to Z = 1. The calculations were 
mainly implemented for the time evolutions in the range of t — 0.5 ~ 10 corresponding 
to the real time range 1.5 ~ 30 million years and the real range in thickness is 0.5 km ~ 
10 km which is the one of main interest in the petroleum industry. In addition, the 
timescale can be chosen in such a way that t = 0.5 ~ 10 corresponding to the real time in 
the order of 15 days to 20 years with a real thickness from 5 ~ 1500 m in civil engineering. 
Numerical results are briefly presented and explained below. The comparison with the 
asymptotic solutions for equilibrium state will be made in the next section. 

Figure 1 shows the poroelastic compaction profile of porosity <fi versus the rescaled 
height Z at different times t = 1,2, 3, 5, 8. The value of A = 100 has been used in the 
calculations. We can see that porosity decreases quite dramatically at the top, and profile 
is nearly exponential versus the rescaled depth 1 — Z. 

Figure 2 provides the viscous compaction profile of porosity versus the rescaled height. 
All the other parameters are the same. The only difference from that of Figure 1 is that 
the compaction relation is now viscous. Comparing with the profile in Figure 1, it is 
clearly seen that porosity changes less slowly than that in the poroelastic case. The 
profile now is more or less parabolic. Although these two figures are quite different in the 
top region, there are still some similarity in the lower region, where the porosity decrease 
very slowly due to the fact that permeability k((p) = (0/0o) m is getting virtually very 
small as < 0o an d m — 8, which will in turn constrain the density- driven flow through 
the porous media, and thus consequently slow down the compaction process. 

To understand these phenomena and to verify these numerical results, it would be 
very helpful if we can find some analytical solutions to be compared with. However, it 
is very difficulty to get general solutions for poroelastic compaction equations (24) and 
(25) or viscous compaction equations (28) and (29) because these equations are nonlinear 
with a moving boundary h(t). Nevertheless, it is still possible and very helpful to find 
out the equilibrium state and compare with the full numerical solutions. 

4.3 Equilibrium State 

To find out the solutions for the equilibrium state, we must solve a nonlinear or a pair of 
nonlinear ordinary differential equations whose solution can usually implicitly be written 
in the quadrature form. In order to plot out and see the insight of the mechanism, we also 
need to solve these ODEs numerically although the solution procedure is straightforward. 
However, it is practical to get the asymptotic solutions in the explicit form in the following 
cases. 

4.3.1 Poroelastic Compaction 

For the poroelastic compaction, the equations for equilibrium state become 




(35) 
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Substituting the expression for p and integrating the above equation once together with 
the top boundary condition (26) gives 

VI -*<H^ -U = (*- *)(!-*.), (36) 

where we have assumed that rh(t) = 1 and h = const. The solution of this equation can 
be written in a quadrature although it is nonlinear. 

Since A = 0.01 —1000, we can expect that two distinguished limits A — > and A — > oo 
will have very different features. For A — > 0, we have 

h — m, ~ 0o, (37) 

which means that porosity does not change and no compaction occur. This corresponds 
to the case of very fast sedimentation or the density difference Ap = p s — p% — > 0. On 
the other hand, as A — > oo, we have 

l±g - 1] » 0, (38) 
a<p oz 

its solution with the top boundary condition can be straightforwardly written as 

= 0oe" a(/l - 2) , (39) 

which is essentially the Athy's profile derived from real field data in sedimentary basins. 
Clearly, if a — > (very slow consolidation), fa <f> , which means that porosity changes 
also very slow. If a — > oo (very quick consolidation), (p — > for h — z > 1/a, which 
implies that compaction proceeds so fast that the porosity is virtually zero everywhere 
except in a thin boundary region at the top. The thickness of the top boundary layer is 
approximately 1/a, which is usually 0(1). However, the solution (39) also satisfies the 
bottom boundary condition |^ — a<p = at z = 0, which means that this solution is a 
uniformly valid solution for steady state. 



4.3.2 Viscous Compaction 

For the viscous compaction, the equilibrium state is governed by 

"= A K r[ l- (i -«- (4o) 

The integration of the first equation together with the top boundary condition leads to 

d (m-/Q(l-fo) 
P dz [ 1-0 J ' 1 } 

and 

= H ± n{m _ _ M ^ _ (1 _ (42) 
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whose general solution can also be written in a quadrature. However, two distinguished 
limits are more interesting. Clearly, if A — > 0, we have 

h = rh, = 0o, (43) 

which is the case of no compaction as discussed in the case of poroelastic compaction. 
Meanwhile, if A — > oo, we have 

(rh - h)(l - W^^) - (1 - 0) = 0, (44) 
which can be rewritten as 

(m-A)(l-0 o )<-^ = O, ^ = Y^- (45) 

By using ip" = ipdip' / dip and integrating from h to z, we have 

(*-*>(l-*0 W . = h ,.» (46) 
2 ip l-0o 

Further integration leads to 



ijerf erf 



2(1 -¥(k- (47) 



1-00 \ 7r(m - /i) 

The comparison of poroelastic solution (39) and viscous solution (47) with the numerical 
results is shown in Figure 3 in the top region where the compaction profile is nearly at 
equilibrium state for A = 1000 and t = 10. The clearly agreement verifies the numerical 
method and the asymptotic solution procedure. 



5 Discussions 

Conventional studies of compaction in porous media have focused on the separate features 
of poroelastic and viscous compaction. The novelty of this paper is to compare and find 
out distinguished features of these two different compaction styles. 

Based on the pseudo-steady state approximations, the model equations of compaction 
can be simply written in dimensionless form as a mass conservation and Darcy's law. A 
constitutive compaction relation is needed to complete this model. In the case of poro- 
elastic compaction, we use an Athy-type relation p = p(0); while in the case of viscous 
compaction due to pressure solution creep only, we choose p = — These two different 
relations result in two quite different behaviours of porosity evolution. In the simpler 
poro-elastic case, we have a single non-linear diffusion equation for porosity 0. 

The analysis showed that the limit A — > (very slow compaction) can be simply 
analysed by means of a boundary layer analysis at the sediment base. The more in- 
teresting mathematical case is when A >> 1 (fast compaction). For sufficiently small 
times, the porosity profile is exponential with depth, corresponding to an equilibrium 
(very long time) profile. However, because of the large exponent m in the permeability 
law k = (0/0 o ) m , we find that even if A >> 1, the product \k may become small at 
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sufficiently large depths. In this case, the porosity profile consists of an upper part near 
the surface where Xk » 1 and the equilibrium is attained, and a lower part where 
AA; << 1, and the porosity is higher than equilibrium which appears to correspond accu- 
rately to numerical computations. For the case of viscous compaction, porosity reduction 
occurs throughout the basin, and the basic equilibrium solution which applies near the 
surface is a near parabolic profile of porosity. The differences in these two profiles are 
very distinguished. 

From the solution (39) for poroelastic compaction at equilibrium state, we see that 
(f) <C 0o when a(h — z) — 0(1) or (h — z) = 0(l/a), that is to say, the solution is 
significant in a region shallower than 

U p « i, (48) 

which corresponds to a depth of 1000 m when d 1000 m and a = 1.0. On the other 

hand, the viscous solution (47) only becomes significant when \Jii{m — h)/2(l — O ) = 
0(1), or in the region of depths h — z greater than 



n 



V 



V T~ h l (49) 



which is equivalent to a depth of 970 m with values of 0o — 0.5, rh — h = 0.3 and 
d = 1000 m. Therefore, we can generally anticipate that the poroelastic compaction is 
dominant in the shallow region from the surface to a depth of 1 km. At depths greater 
than 1 km, the pressure is high enough, pressure solution mechanism becomes significant 
and thus compaction is essential viscous. Naturally, there exists a region of depths near 
lkm where both mechanism becomes important, and an obvious extension is to include 
both models in a more realistic model. From the poroelastic constitutive relation (4) 
and viscous relation (8), we can formulate a generalised viscous-poroelastic compaction 
model of Maxwell type 

(50) 

Subsequently, we would expect a visco-poroelastic porous medium and thus some care 
is needed to ensure the resulting model involving material derivatives is frame invariant. 
Fortunately, this frame invariance is alway true in the present 1-D formulation. Incorpo- 
ration of these extension and other processes such as convection and 3-D density-driven 
flow will form the substance of future work. 

Acknowledgements. The author wishes to thank the anonymous referees for their 
very helpful comments and very instructive suggestions. 



References 

1 I. Lerche, Basin Analysis: Quantitative Methods ( Academic Press, San Diego, California, 1990). 

2 R. E. Gibson, G. L. England and M. J. L. Hussey, The Theory of One-dimensional Consolidation of Saturated 

Clays, I. Finite Non-linear Consolidation of Thin Homogeneous Layers, Can. Geotech. J., 17(1967)261- 
273. 



10 




0.1 0.2 0.3 0.4 0.5 

Porosity 
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Figure 2: Viscous compaction profile of porosity versus the rescaled height Z. All the other parameters 
are the same as in Figure 1. A viscous compaction relation between effective pressure and velocity is 
used. The profile now is nearly parabolic. 
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Figure 3: Comparison of asymptotic solutions (39) and (47) (dashed curves) with numerical results (solid 
curves) in the top region where the profile is nearly at equilibrium state. The agreement is clearly shown. 
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